A precise approximation for directed percolation in d = 1 + 1 
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Abstract 

We introduce an approximation specific to a continuous 
model for directed percolation, which is strictly equivalent to 
1 + 1 dimensional directed bond percolation. We find that 
the critical exponent associated to the order parameter (per- 
colation probability) is (3 = | (l - -±^J = 0.276393202..., in 

remarkable agreement with the best current numerical esti- 
mate (3 = 0.276486(8). 

PACS numbers: 05.70.Ln, 64.60.Ht, 64.60.Ak 

Directed percolation (DP) Jl|-Q is a useful paradigm 
for dynamical phase transitions between an ac- 
tive/spreading phase and an extinct /absorbing phase. 
Models in the DP class of universality are involved in 
the description of catalytic reactions ||, surface dynam- 
ics ||, porous systems g, granular media [Q, epidemics, 
Calcium dynamics in cells [Q, developed turbulence and 
coupled maps [^)... Recently, H. Hinrichsen summarized 
the large scope of possible physical applications of DP 
H, which led P. Grassberger to conjecture that the DP 
universality class should describe any continuous phase 
transition from a fluctuating active phase into a single ab- 
sorbing phase, in the absence of quenched disorder and 
special symmetries . In a sense, DP plays a similar role 
in the study of dynamical phase transitions as the Ising 
model for continuous equilibrium phase transitions [jjj. 

Despite its ubiquity, DP is maybe the only major sta- 
tistical physics model which has not yet been successfully 
solved in one spatial dimension (+ time), probably due 
to its lack of conformal invariance. 

Let us recall the original model of directed bond per- 
colation in d = 1 + 1, describing the propagation of a 
fluid in a 2c? porous medium. On a square lattice tilted 
at 45°, a fraction of bonds p are chosen at random to be 
active, whereas the remaining bonds stay inactive or bro- 
ken. The "fluid" starts from the top row, and propagates 
downward, only passing through the active bonds. One 
then defines the order parameter n(p,t), which measures 
the average density of occupied sites at row t. nip, t) 
happens to coincide with the probability that at least 
one site at row t is still active (percolation probability). 

In the stationary limit t — > +oo, the order param- 
eter tends to a constant value nip), which is zero be- 
low p c = 0.644700185(5) $3^, and behaves as n(p) ~ 
(p — p c )P , near p c . This defines the universal critical ex- 
ponent = 0.276486(8) 
m{t) — 0), when the i th 
(resp. inactive/empty), 



|0|Jll). Defining ^(t) = 1 (resp. 
site at row t is active/occupied 
rii(t) satisfies the following re- 



cursion relation 

1) = a^iit) +6 i n i+1 (i) - aA^iW^i+itX), (1) 

where cii(t) and bi(t) are independent random variables 
taking the value 1 with probability p (if the corresponding 
ongoing bond is active), and otherwise. 

In relation to self-organized criticality (SOC), it has 
been recognized that directed bond percolation is strictly 
equivalent to a continuous dynamical model (SOCDP), 
involving no external parameter (like p in DP) [|l2|— (l5| . 
On a Id lattice, we define the continuous variables Xi(t) 
as satisfying the recursion relation 

Xi(t+ 1) = mm[max(Xi(t), Zi),Taax(xi+i(t), z'i)] , (2) 

where £j(t) and z[{t) are independent random variables 
uniformly distributed between and 1. It can be easily 
shown that the n^'s for directed bond percolation and the 
Xj's are very simply related: 



n,-(t) = 0(x l {t) -p), 



(3) 



where 6(.) is the usual Heaviside step function. In the 
large time limit, the x^s are distributed according to a 
stationary probability distribution p(x), and 



n{p) 



p{x) dx. 



Now, using Eq. (Q) (or equivalently Eq. 
the stationary limit, 



2p-l 



nip)-- 



2p- 1 



(nin 2 ) 



P2(xi,x 2 ) dxidx 2 



(4) 

we find in 

(5) 
(6) 



where p 2 ixi, x 2 ) is the nearest neighbor correlation func- 
tion of the x,;'s. In mean field (MF) theory one makes 



the approximation {n\n 2 ) 
2p-l 



(ni) , leading to 



riMFip) = 



V 



Pmf(p) 



2(1 -p) 



(7) 



From now, we study the stationary state of DP in terms 
of the continuous model defined by Eq. (||). We first 
notice that 

Xi(t + 1) = minfa;!^), x^it)], with probability (8) 

Pmin = min[xiit),x 2 it)], 
xi(t + 1) = max[xi(t), x 2 (t)], with probability (9) 

Pmax = max[ari(i), a: 2 (*)](! - max[xi(t),x 2 (t)}), 



1 



so that there is a non zero probability that x\(t + 1) = 
x 2 (t + 1) exactly. Hence, the two-point correlation func- 
tion of the Xi's, p 2 (xi,x 2 ), should include a 8(x\ — x 2 ) 
contribution (<5(.) is the Dirac peak distribution). Thus, 
in all generality, we write p 2 {x\, x 2 ) in the following form: 

p 2 {xi,x 2 ) = p 2 (xi,x 2 ) + p(xx)g(xi)6(xx - x 2 ), (10) 

which defines g(p) as the probability that x 2 = p, con- 
ditional to the fact that its neighbor x\ = p. We then 
define f(xi,x 2 ) through the relation 



p 2 {xi,x 2 ) = p(mm(xi,x 2 ))f(xi,x 2 ), 



(11) 



noting that as p(p) diverges near p c (since (3 < 1), 
p(min(xi, x 2 )) > p(max(xi,x 2 )), at least near p c (numer- 
ically p(p) appears to be a strictly decreasing function, 
like in mean field theory). We expect that f(xi,x 2 ) is a 
smooth function of order unity. Indeed, contrary to the 
MF approach (where one assumes that (nin 2 ) ~ [n(p)} 2 ), 
correlation functions all behave as n(p) near p c . Indeed, 
as p c > 1/2, Eq. (JsJ) implies that (nin 2 ) ~ n(p) ~ 
{p-PcY, such that p 2 (p,p) ~ p{p) ~ (p - p c )~ (w,) 
(instead of p(p) 2 ~ (p — Pc)~ 2 ' 1_,3 \ predicted by MF). 
A natural guess for f(xi,x 2 ) is provided by the general 
statement that although MF is inept at describing cor- 
relation functions near p c , it still leads to reasonably ac- 
curate amplitude ratios between them, for all values of p 
(at least for short range correlation functions like (nin 2 )). 
Hence, this prompts the introduction of the key approx- 
imation 



f(xi,X 2 )^ fMF(xi,X 2 ), 

_ Pmf(xi)pmf(x 2 ) 

p M F(min( x i, x 2)) ' 
= Pmf(^bx(x 1) x 2 )). 

In the following, we make the more general ansatz 

P 2 (xi,x 2 ) = p(min(a;i,a;2))/(max(a;i,a;2)), 



(12) 

(13) 
(14) 

(15) 



where f(p) is not necessarily equal to Pmf(p)- pip), fip) 
and g(p) are not independent functions as they are re- 
lated together by Eq. (||), and by the probability conser- 
vation constraint, 



pip) 



p 2 {x,p) dx. 



(16) 



From Eq. (|5|) and Eq. (jig), and after straightforward 
calculations, we obtain the two relations 



n(p) 

g(p) - 



exp 
i 



p 



1 2f(x) - pmf{x) 
n M F(x) - g(x) 



f( x )dx + f(p)^- = l. 

Pip) 



dx 



(17) 
(18) 



From Eq. ( |i"7| ) and Eq. (jll), one can obtain a first order 
differential equation for F[p) — f(p) — Pmf{p), involving 



only p and g(p). This equation can be shown to have 
only F = as a global solution satisfying the boundary 
conditions and the physical constraints. Hence, and in 
complete accordance with the physical argument given 
in Eq. (O), we find 



fip) = PMFip)- 



(19) 



Quite remarkably, for this precise form for f(p), Eq. dl§| ) 
is now satisfied for any choice of g(p), so that we are left 
with Eq. (U7]) as the only non trivial relation between 
n{p) and gyp)- 

Now, as n(p) vanishes at p Cl we expect the function 
involved in the integral of Eq. (]l7]) to develop a single 



pole at pc, of residue /3, so that n{p) 
p c - This leads to 



ip - Pcf 



liPc 



nMFKPc) = 2" 



g'(Pc) 
PMFiPc) 



1 



(20) 
(21) 



Note that Eq. ([20|) is in fact an exact identity, which does 
not rely on the present approximation on p 2 (x,y). 

In order to achieve our goal of computing n(p), we 
need a further relation for g(p). This can be obtained by 
writing the exact stationary equation for g(p). We first 
define (•••) as the probability of having three consecutive 
sites with xi = p, divided by p(p) (that is, conditional to 
having one site with X{ = p). In the same manner, we 
define (• • J>) (resp. (• • p)) as the probability of having 
two consecutive sites with x\ = x 2 — p, and X3 > p (resp. 
X3 < p), divided by p(p). Finally, in the stationary limit, 
we obtain the exact relation 



g{p) = {2p-p 2 ) 2 {...) 



(22) 



+2{2p-p 2 )(p(mmp) +p(l -p) (..p)) 
+p 2 (p*p) +2p 2 (l -p)(p»p) +P 2 il-p) 2 (p»p) 



+P 2 (»p») 



For instance, the first term (2p—p 2 ) 2 (•••) represents the 
fact that a configuration • • (x\(t + 1) = x 2 (t + 1) = p) 
at time t + 1 can arise from a configuration • • • {x\ (t) — 
x 2 {t) = x 3 (t) — p) at time t, provided that x\ and x 2 are 
preserved by the transformation of Eq. (||) . This happens 
with probability 

(Pmin+Pmax) 2 = ip + p(l - p) f = (2p - p 2 f \ (23) 

hence the coefficient in Eq. d23) (pmm and p m ax have been 
defined in Eq. (|) and Eq. If). 

Eq. (|22] ) relates g(p) to three-point correlation func- 
tions, and cannot be exploited unless an additional ap- 
proximation is introduced. We will factor these three- 
point correlation functions into products of two-point 
correlation functions, according to the usual mean field 



2 



scheme. Introducing p + as the probability that x-i > p, 
conditional to the fact that x\ = p, we obtain 



Jp h(x,p)dx 



[i - g(p)]p+= 



V 



pip) 

Pmf(x) dx = 5 . 

P 



(24) 
(25) 



where Eq. (|2J) is an exact identity. We give below a few 
examples of three-point correlation functions computed 
according to this MF factorization scheme: 



(•••>=. 9(P) 2 , 

(• •P) = fl(p)[l - ff(p)]p+ = ff(p) 
= [1 - ,g(p)] 2 p + (l -p+), 

g{p) 



(i-pf 



2p-l 



P 



{l-pf 



P 



(26) 
(27) 
(28) 
(29) 



Inserting the MF form for the three-point correlation 
functions into Eq. (|22|), we finally obtain a closed equa- 
tion for g(p), 

g( P ) = (l-p + pg( P )) 2 + (l-p) 2 g(p)(l-g( P )), (30) 
which can be readily solved, leading to 

{l-pf 



g{p) = 



2p-l 



(31) 



g(p)=(1-p)Wl) 
g(p) from NS 




FIG. 1. We plot (• • •) (two top lines) and (• • p) as ob- 
tained from numerical simulations (full lines), and as given 
by Eq. ( p6i|) and Eq. (E?]) (dashed lines), where the numerical 
value of g(p) has been inserted in these expressions. Note 
that (• • •) ~ g(p) 2 , especially near p c . These functions 
all vanish as (1 — p) 4 near p = 1, as predicted by Eq. ( |2^ ) 
and Eq. @. Insert: comparison between the numerical 
g(p) and the present theory. In all figures of this letter, we 
have simulated a system of N = 300000 sites, averaged over 
100 samples. Physical quantities in the stationary state have 
been estimated by averaging them between t = 300000 and 
t = 310000. 



p c and f3 can now be calculated by expressing the con- 
ditions of Eq. © and Eq. @. We obtain 



Pc 



g{pc) 



VE- l 



0.618033989. 



(3 = - ( 1 - -4= ) = 0.276393202. 

2 V vV 



(32) 
(33) 



where r is the golden mean. p c is only in fair agreement 
with the best numerical estimate p c = 0.644700185(5) 
]l0| , p| , although this represents a dramatic improvement 
when compared to the mean field value p c = 1/2. Note 
that the exact identity Eq. ( |20| ) implies that getting 
p c > 1/2 necessitates the introduction of a non trivial 
function g(p), which is zero in MF. The (3 exponent is in 
extraordinary agreement (relative accuracy of 0.034%) 
with the best available numerical value (3 = 0.276486(8) 

US 




0.7 0.8 0.9 1.0 
P 

FIG. 2. We respectively plot n(p) as given by mean field 
theory (dashed line), the present theory (dotted line), and 
numerical simulations (full line). Note that the three curves 
coincide near p — 1, as mean field theory becomes exact in 
this limit. A fit of the numerical data to the functional form 
of Eq. (|^), where r becomes a fitting parameter, cannot be 
distinguished from the actual data. 

The fact that the relation (• • •) = g(p) 2 seems to 
be exactly satisfied numerically at p c could explain this 
agreement, which also implies that the MF factorization 
of Eq. (^2|) is quantitatively correct near p c (the three- 
point correlation functions appearing in Eq. (|22| ) are not 
independent and are related to (• • •) by various sum 
rules). This is illustrated in Fig. 1, where the exact nu- 
merical g(p) = (• •), (•••), and (• • p) are plotted with 
their theoretical counterpart. 

Now, f(p) and g(p) being known, the percolation prob- 
ability can be easily computed by using Eq. 



n{p) = p 



(x - r)(2 +t-x) 



/3 



(34) 



{x-1 + t)(x + \+t) 



1 + T 



1-/3 



3 



In Fig. 2, we compare this result with the numerically 
extrapolated stationary percolation probability, and to 
the MF result of Eq. (g). 

Finally, in Fig. 3, in order to test the validity of our 
basic approximation Eq. (fl4|), we plot f(xi,X2) (defined 
in Eq. ([tl])) as a function of max(xi, a^)- We find that 
this scatter plot is reasonably aligned around an effec- 
tive curve, and that /mf(iiiI2) = PMF(m&x-(xi, X2)) 
appears to be a lower bound for the actual f{x\, £2). 
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function of max(xi, 22) (scatter plot). This is com- 
pared with the approximation central to this paper: 
f(xi,X2) ~ /)Mp(max(xi,a;2)) (thick line). We observe that 
the theoretical expression seems to be a lower bound for the 
actual f(xi,X2), and that the dispersion (due to the explicit 
dependence on min(xi, 2:2)) is weak enough, so that the scat- 
ter plot tends to align around an effective curve. 

In conclusion, we have introduced a new approxi- 
mation for a continuous model equivalent to directed 
bond percolation. In this language, this approximation 
amounts to properly modelizing the correlation function 
PiipxiV?) 1 relating the properties of directed bond per- 
colation for two different percolation parameters p\ and 
P2- By assuming that amplitude ratios are correctly de- 
scribed by mean field theory, we end up with a precise 
description of the percolation probability. In particular, 
we find an exponent f3 in remarkable agreement with the 
best available numerical simulations. 

It would be interesting to exploit the present approach 
in order to describe the dynamical properties of DP. This 
study is currently in progress. 

This approach could also prove useful in tackling the 
notably difficult problem of parity conserving branching 
annihilating walks || . This universality class is exempli- 
fied by the reaction-diffusion model of diffusing particles 
A, involving annihilation (A + A — > 0) and branching 
(A — > A + A + A) processes. This problem has so far 
eluded all manner of theoretical approaches in d = 1 + 1. 



[1] J. Marro and R. Dickman, Nonequilibrium phase tran- 
sitions in lattice models, Cambridge University Press, 
Cambridge (1999). 

[2] W. Kinzel, in Percolation Structures and Processes, Ed. 
G. Deutscher, R. Zallen, and J. Adler, Ann. Isr. Phys. 
Soc. 5, 425 (Adam Hilgor, Bristol, 1983). 

[3] H. Hinrichsen, Adv. Phys. 49, 815 (2000). 

[4] H. Hinrichsen, Braz. J. Phys. 30, 69 (2000). 

[5] R. M. Ziff, E. Gulari, and Y. Barshad, Phys. Rev. Lett. 
56, 2553 (1986). 

[6] L.H. Tang and H. Leschhorn, Phys. Rev. A 45, R8309 
(1992). 

[7] H. Hinrichsen, A. Jimenez-Dalmaroni, Y. Rozov, and E. 

Domany, J. Stat. Phys. 98, 1149 (2000). 
[8] Y. Pomeau, Physica D 23, 3 (1986). 
[9] P. Grassberger, Directed percolation: results and open 
problems, preprint WUB 96-2 (1996), unpublished. 
[10] I. Jensen, Phys. Rev. Lett. 77, 4988 (1996). 
[11] I. Jensen, J. Phys. A 32, 5233 (1999). 
[12] A. Hansen and S. Roux, J. Phys. A 20, L873 (1987). 
[13] P. Grassberger and Y.-C. Zhang, Physica A 224, 169 
(1996). 

[14] S. Maslov and Y.-C. Zhang, Physica A 223, 1 (1996). 
[15] R. Dickman, M. A. Muiioz, A. Vespignani, and S. Zap- 
peri, Braz. J. Phys. 30, 27 (2000). 



4 



